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Abstract. We provide an explicit method for solving general markovian master equations for 
quadratic bosonic Hamiltonians with linear bath operators. As an example we consider a one- 
dimensional quantum harmonic oscillator chain coupled to thermal reservoirs at both ends of the 
chain. We derive an analytic solution of the Redfield master equation for homogeneous harmonic 
chain and recover classical results, namely, vanishing temperature gradient and constant heat current 
1-^ in the thermodynamic limit. In the case of the disordered gapped chains we observe universal heat 

current scaling independent of the bath spectral function, the system-bath coupling strength, and the 
boundary conditions. 
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g 1. INTRODUCTION 

I 

pH The question of the microscopic origin of well established Fourier's law of heat conduc- 

O tion has remained unsolved for many decades [1, 2, 3, 4, 5, 6]. The main issue is to find 

i_Zj the necessary and sufficient conditions for normal transport. Therefore, it is important 

to understand basic, simple models, that are expected to embrace important features of 
complicated, more realistic models. It is believed that for solid systems disorder and 
0^ anharmonicity play a key role in determining the transport properties. Focusing on the 

^^ quantum of quantum mechanical systems, this encouraged the development of different 

pq formalisms for treating large non-equilibrium quantum systems, such as the Keldysh 

^ technique [7], Landauer-Butikker formalism [8], quantum Langevin equation [9], quan- 

ta tization in the Fock space of operators [10, 11, 12, 13] giving insight into the region far 

CN from equilibrium. 

The simplest, but nevertheless non-trivial model in the context of heat transport 
. ^ properties, is the coupled one-dimensional harmonic oscillator chain. Many different 

^ methods on modeling the heat reservoirs have been presented. Much effort has been 

given to solve the classical many-body Langevin equation, or the corresponding master 
(Liouville) equation. For the homogeneous harmonic oscillator chain size independent 
heat flux and vanishing temperature gradient have been proven [14]. On the other hand, 
a model dependent scaling of the heat flux has been found for the disordered chain [15]. 
Dhar and Shastry have derived a generalized quantum Langevin equation by using the 
Ford-Kac-Mazur formalism [16], and have got similar results as in the classical case. 
Further, it has been argued in detail that the heat and particle transport properties of the 
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model depend on the spectral properties of the heat bath [6] and boundary conditions of 
the model, in the case of a gapless harmonic oscillator chain, i.e. with vanishing on-site 
potential. In Ref. [17] ordered and disordered one-dimensional harmonic chains within 
the quantum mechanical Langevin equation have been analyzed. No clear statement on 
scaling of the heat current in the disordered chain has been given. The connection to the 
entanglement has been observed as well. 

In this paper we consider a different, effective approach of open quantum systems 
[18] in the many-body context, which in the past lead us to interesting and fruitful 
results in the realm of one-dimensional fermionic models [10, 11] and can be used 
to investigate heat transport properties via the Redfield master equation [19, 20]. It 
has been shown [13] that the Lindblad master equation, a special case of the Redfield 
master equation, can be solved exactly for quadratic bosonic models. We shall see that 
similar formalism can be used to get the non-equilibrium stationary state (NESS) of the 
Redfield master equation for a general quadratic boson Hamiltonian with linear bath 
operators. We further apply the method to one-dimensional harmonic oscillator chain, 
which serves as a well known model and can be used to test our approach and possibly 
give an alternative description of heat baths. 

This paper is organized as follows. In section 2 we review the main concepts of the 
quantization in the Fock space of operators for bosonic systems [13]. We show that 
this method can be used to solve the Redfield master equation and derive the Lyapunov 
equation for the covariance matrix, which determines the NESS. In section 3 we consider 
the one-dimensional harmonic oscillator chain model in the context of the Redfield 
master equation. First, in subsection 3.1 we describe a general method of solution of 
the open harmonic chain and in subsections 3.2, 3.3 provide two explicit solutions for 
the homogeneous chain, namely a compact solution for the case of symmetric (left/right) 
system-bath coupling and ohmic bath correlation function and a recursive perturbative 
(in the coupling) solution, again for symmetric coupling strengths and a general bath 
correlation function. In both cases simple expressions for the heat conductivity are 
obtained. Further, in subsection 3.4 we introduce disorder and find power law scaling of 
heat current, that is independent of the fine reservoir properties, i.e. the spectral function 
of the bath, and the coupling strength. In section 4 we summarize and conclude. 

This paper mainly presents original research results which have not been published 
before, however we review some material from [13] and [10, 11] in order to make the 
manuscript more smoothly readable and self-contained. 

2. EXACT SOLUTION OF THE REDFIELD EQUATION FOR 
QUADRATIC BOSON SYSTEM 

Our aim in this section is to extend operator Fock space quantization for bosonic open 
quantum system. We focus on the Redfield master equation since the Lindblad case can 
be considered as a limit of equal-time bath correlation functions similar as in Ref. [11]. 
We shall provide an exact solution of the Liouville master equation in the Redfield form 

^p(0 = ^pW, (1) 

^p = i[p,//]+#p, (2) 
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^P =LMv/o°°d?rv,^(T)[X^(-T)p,Xv]+h.C. (3) 

The Liouvillean (2) includes two parts. The first part is the generator of the unitary 
evolution and is determined by the Hamiltonian H. The second part, the dissipator, 
introduces effectively the influence of infinite large environments (reservoirs or baths) 
with the correlation functions Yy^pi{t). We assume a general quadratic Hamiltonian and 
linear coupling operators 

// = p.Pp + ^.Q^ + p.R^ + ^.Rp, (4) 

Xy ^ Xy ■ q -\-Xy ■ P, 

where P, Q, R are real, symmetric matrices, Pj.qj are momentum and coordinate opera- 
tors, respectively, and x^ -.Xy ■ are real numbers (coupling amplitudes), which determine 
the coupling operators Xy. Greek subscript letters are used to denote different reservoirs 
and latin subscript letters denote the position in the system (site index j = 1 , 2, . . . n) 
of size n. We shall use this notation throughout the paper, as well as bold latin letters 
for matrices (A) and the hat for super-operators (a). In the equation (4) we introduced 
the following notation of a column vector x = (jci,X25 • • -Xn)^, where x can represent a 
number, an operator or a super-operator, which is evident from the context, as well as 
a dot product that assigns only transposition and not the complex conjugate of the left 
operand x-y = L/=i -^jj; • 

The problem could equally well be formulated in the second quantization picture 
using the creation and annihilation operators a], a:. It is obvious that the two representa- 
tions are equivalent, however, using the Hermitian operators results in a real Lyapunov 
equation for the stationary correlation matrix. This simplifies the algebra and the numer- 
ical calculations. 



2.1. Coordinate-momentum picture super-operators 

We conform to Ref. [13], where the solution of the Lindblad master equation has 
been found exploiting the algebra in a pair of dual vector spaces J^ and J(^', where ^ 
contains trace class operators (density matrices) and J^' contains unbounded operators 
(physical observables) ^ . We will adopt the Dirac notation and write an element of J^' 
as ket \p) and an element of ^ as bra (A|, so that their contractions give the expectation 
value of A for the state p 

(A|p)=tr{Ap}. (5) 

The basic tool used in this section are left and right multiplication maps over ^ defined 
as 

b'^\p) = \bp), b^=\pb). (6) 



For details see Ref. [13]. 
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From the cyclicity of the trace and the definition of the inner product (5) we deduce the 
action of the maps (6) on the elements of the adjoint space 

{A\b^ = {Ab\, {A\b^ = {bA\. (7) 

Further we define An maps Cyj.Cy for j = 1 , 2, . . . , n and v = 0, 1 

^Oj = p). £'oj = Kqyqj). (8) 

satisfying almost canonical commutation relations 

[c^j, cv,^] = [c^j, Cy^i^] = 0, [c^j, Cy^i^] = 5j^k5^,v (9) 

Note that c' ■ 7^ c!^ ,• Another important property is that the maps Cy left- annihilate 
the identity operator (l|c' = 1. We shall now formulate the Liouvillean as a quadratic 
form in the maps Cy /, c' using the inverse relations 

Pj=(^0,j, q) = Cx^j-lCQ., (10) 

Let us first rewrite the contribution of the Hamiltonian to the Liouvillean 

i(^'^ - H^) = i{p^ ■ Vp^ + q^ ■ Q|R + q^ ■ Rf^ + p^ ■ Rq^ (11) 

= 2{c[ ■ Pco - c[) ■ Qci + c[ ■ Re, - c(, ■ Rcq) + i{-c[ ■ Pc[ + c^, ■ Qc'q) . 

Above, anti-Hermitian part of the Liouvillean (11) left- annihilates the identity operator 
and is quadratic in the maps Cvj.c'y ■. We shall shortly see that these characteristics are 
satisfied for the second, dissipative part of the Liouvillean (3) as well. The dissipator can 
be represented in the following, useful form 

^ = lT! t <,. [d<j(-') (r5,M(^)4'r +r5,;(T)4;) . (12) 

^l,V l,m j,k=l JO ^ 

We introduced the Heisenberg propagator 

Xy = {xl,xl), fl{t)=xl-exp{-iadHt), aG{p,q}, (13) 

and the fundamental dissipators 

'^J:k\P) = \[^j^Ph]), ^f,\p) = \[bkp,aj]), j,k=\,...n, (14) 

where aj = pj if a. = p or Uj = ^^ if a = q and bj = p j if b = p or bj = Qj if b = q. 
Expressing the relations (10) we obtain after some algebra a compact form of the 
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fundamental dissipators 

-^fS = '<,kici,j - K,j)^ ^f^ = -ic[,/ij, (15) 

Let us now rearrange and group the terms in the dissipator (12) according to their 
meaning. All information about the baths can be encoded in the bath matrices 

M'"'' = £ rdTte®/J^(-T))r5,^(T), /,mG{p,q}, (16) 

which simplifies the equation (12) for the dissipator 

By (•)* we denote the complex-conjugation. Plugging the relations (15) in the equation 
(17) we obtain an apprehensible, quadratic form of the Redfield dissipator 

^ =2{-c^-Mf''^CQ-2c'Q-Mf%+c[-Mf'^CQ + c!i-Mf%) (18) 



We shall use the following notation Mr'*" = Re(M''''') and Mf'*' = Im(M^''') for real and 
imaginary part of the matrices M^'^, a,b G {p,q}. Inserting now the forms (11) and (18) 
in the definition (2) and defining the super- operator vector h = (co,Ci,Cq,Cj)^ we arrive 
at the compact form of the complete Liouvillean 



^ = b-Sb- 


-Sol, 


s ( « 


-X\ 
Y J' 


■ R + M^i'P, 
-P-Mf'P, 


-R-Mf''i 



(19) 



Y-lf 2iQ + Mi'i + (Mi'i)T, _(Mq'P)*-(MP'i)T 
2 1^ -(Mi'P)t - (MP'1), -2iP+ (MP'P)* + (MP'P)'^ 
5o = 2Im(Tr(MP'i-Mi'P)). 

We name the part of the matrix X arising from the dissipator the coupling matrix and 
the matrix Y the driving matrix. The left stationary state of the Liouville equation (1) is 
simply the identity operator, (l|=if = 0, whereas the right stationary state determines 
the density matrix Pness of the non-equilibrium stationary state, =SfpNESS = 0- We 
now employ the main theorem of Ref. [13]. It states, that the right stationary state 
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of a Liouvillean, which may be represented as a quadratic form in the creation (c' ) 
and annihilation (cj) super-operators, that satisfy canonical commutation relations, is 
determined by the covariance matrix (two point correlation function) 

Z,-, = (l|cVfc|NESS). (20) 

This correlation function is calculated from the continuous Lyapunov equation 

X'^Z + ZX = Y. (21) 

The NESS is a Gaussian state, hence all expectation values can be calculated from the 
two point correlation function (20) using the Wick theorem. In our case the correlation 
matrix Z gives us the momentum-coordinate correlations 

/ ZP'P (Z^'P)^ \ ab 

Z=( ^q.p' Zi'1 J' 2;>=tr{«APNESs}, a,bG{p,q}. (22) 

The correlation matrix Z is made real by adding the matrix Zq := ^(j'^'^ln- The map 
Z -> Z -(- Zo implies Y -)■ Y - X^Zq - ZqX, therefore we have 

/ ZP,P^ (Z'l'P)^ \ ~ i 

'^~\ zq,p' zq,q )' Z ' =Z ' +2^n, (23) 

-, _ 1 / M?''! + {Ury, -M?'P - (MP''^)T 
^ " 2 V -(MD^-MP'^ MP'P + (MP'P)t 

The coupling and the driving matrices are connected to the real and imaginary part of 
the bath matrices (16), respectively. The Liouville super- operator in Redfield form does 
not necessary conserve positivity, therefore one should check if the obtained covariance 
matrix determines a bona fide density matrix. This is done by verifying the positivity 
of the covariance matrix Cjj^ = ti:{bjbkp^ESs}^ which is simply an unusual form of the 
Heisenberg uncertainty relation expressed for the covariance matrix [21]. We introduced 
an operator vector b = {p,q). 



3. HARMONIC OSCILLATOR CHAIN 

In this section we use the derived theory to solve the example of harmonic oscillator 
chain coupled to heat baths at the ends of the chain. The classical and quantum version 
of this problem have been studied in the literature [6, 17] using the generalized quantum 
Langevin equation, derived in Ref. [16] utilizing the Ford-Kac-Mazur formalism. A 
boundary condition dependent scaling of the heat flux for the gapless disordered chain 
has been observed. In Ref. [17] the connection between heat current and entanglement 
has been analyzed. Here we propose a different approach to introduce the effect of the 
heat baths. We use the slightly modified Redfield dissipator, namely, we extend the lower 
integration limit in the equation dissipator (12) to — oo, which is the same as neglecting 
the principal value of the integral. This further simplification is necessary for the Gibbs 
state to be the stationary state, when all baths have equal temperature. It is not clear 
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under which conditions this approximation is valid, however it is necessary to obtain 
physically meaningful heat reservoirs (see Ref. [11] for more discussion on this issue). 
The Hamiltonian of a harmonic chain with nearest neighbor interaction may be 
expressed as a sum of local Hamiltonians and some boundary Hamiltonian 



n-l 



(24) 



The local Hamiltonians Hj contain the contribution from the kinetic energy, on-site 
potential, and the interaction energy, whereas the boundary Hamiltonian consists of the 
kinetic and the potential part on the ends of the chain and the boundary condition term 



Hj = \ {pj + PJ+i + COJq^j + O^Wj+i) + 2 



1j Ij+l 



Hb = \{p\ + pI + (0lq\ + 0)2^9 + ^^q\ + ^jl 
Employing a compact matrix notation we get 

H = jP-p + jq-Qq 



/^i±^ + «2 



Q 



ki 



mi 



^/tnimi ■ 



^m\m2 ' 
m2 + "^2 ' 







^n- 



V 



0, 



^«— 1 



yjnin-imn' 



^mn-\m„ 



We consider a natural choice of the coupling operators 



(25) 



(26) 



Xl = \/£Lqi, ^R = A/CRfJ'n, Cr 



m„ 



£r 



nin 



and two baths with possible different temperatures and general spectral functions 

T^Uf.^-^ sign(«)|«r 



'exp(©/3j)-l 



(27) 



(28) 



Three different types of the heat bath according to the parameter v are possible: a) sub- 
ohmic bath (v < 1), b) ohmic bath (v = 1), and c) super-ohmic bath (v > 1). We use 
the values v = ^ and v = 2 for the sub-ohmic and the super-ohmic case, respectively. 
The main goal of this section is to examine the behavior of heat current defined by the 
continuity equation in the bulk ^ 



J 



j-^J 



= ^tr(//jPNESs) = tr{i[//,//;]pNESs} = Jj-lJ + Jj+lJ, 

tr{i[//;_i,//,]pNESs} = -^tr 



(29) 



Pj<ij-2-Pj<lj-l 



nij^i 



PNESS 



The reasoning leading to equation (29) for the heat current is the same as in [1 1]. 
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where Jj-ij denotes the current from site j — 1 to site j. From the antisymmetry con- 
straint for the correlation matrix Z^p, which we shall derive shortly, follow vanishing 
expectation values {pjQj + qjPj) = 0. Therefore the equation for the heat current sim- 
plifies to 






(30) 



In the remaining part of this paper we investigate the problem of heat transport through 
the harmonic oscillator chain in the framework of third quantization for bosonic Redfield 
master equation. First, we recover some results for the homogeneous chains, then we 
focus on the scaling law of the heat current in the disordered chain. 



3.1. Homogeneous chain 



In this section we study the homogeneous chain, where all coupling constants, masses 
and on-site potentials are equal; kj = A:, my = m, (Oj = (Oq for all j . Thus, only two 
parameters determine the behavior of the system, namely the on site potential (Oq and 

the coupling frequency (Dq - 



. The Hamiltonian is 



H 






„2„2 



f + ^]+LPl4(^7-^,+i)^ + 4(^? + ^^) 



(31) 



Q = (o2 



H = jp-p + jq-Qq, 
( 2 -1, ■■■, \ 

-1, 2, ••., ; 

:, ■••, ■••, -1 
V 0, ■■■, -1, 2 / 



+ Oil^n- 



In this section we restrict ourselves to fixed boundary conditions A:j = k'^ = k = 1, 
whereas in later discussion of the disordered chain we refer also to free boundary 
conditions, k\ = k'„ = 0. The coordinate part of the Hamiltonian is diagonalizable by 
sine-Fourier transformation 



U.. 



^, sm ( ^ 

n+l \ n+i 



Q=u^u^ 



(32) 



Normal coordinate and momentum operators are obtained by the transformations q' = 
q -IJ, p' = p -IJ. The Hamiltonian is then diagonalized via the usual transformation 

to creation and annihilation operators q'- = {a- -^aj)! ^/TXj^ /?'■ = \^fkj{ci- —aj)/\/2. 
We can make this transformations for a general harmonic chain where the columns of 
the matrix U are the right eigenvectors and Xj are the corresponding eigenvalues of 
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the matrix Q. For the homogeneous chain we have U^ = U^ = U and Xj = ^jj- The 
Heisenberg propagator in the normal basis contains a coordinate and a momentum part 

ci'j{t)=ff{t)q'j+ff{t)p) (33) 

J J J 

The time dependent coupling operators in the normal basis are 

^) =\fWi'^^{lSi)^ ^lW =^L-(diag(/'i'(0)^' + diag(/P'(0)/), (34) 
'J =\/1^^in(S)' ^rW =^^-(diag(/i'(0)£' + diag(/P'(Oy)- 



I I OAll I III 

«+l \n+l J 

The problem now is to solve the continuous Lyapunov equation 

XTZ + ZX = Y, (35) 

^, ]' H 0, 



^/„/ „/„/ 



e'^^J-l\-rB,,(^\\ , i^r^R^^Rj- „//e'^R^y-l 



Bath matrices M^i ,M1P are obtained from the equation (16) applying the Fourier 
transformation r^(a)) = j^J°°^dtexp{—icot)rP{t) and the Kubo-Martin-Schwinger 
(KMS) condition r^(-(o) = exp(j8©)r^(©) 

M')'^' = fx^^xi^diag ((l +eA^^^) T'^K^)) + fx^^^ix^diag ((l +e^^i) r^(A;)) , (36) 

In the case of equal temperatures of the baths the Gibbs state should be a stationary state 
of the modified Redfield master equation. Therefore, it is a good exercise and consis- 
tency check to calculate the equilibrium solution (for /3l = j8r = j8) of the Lyapunov 
equation (35) 

Xj(l+eP'A ,, (l+e/5^. 



jyjq'p' ^ I^L ^ ^L^jag f f e^^W-1 pft^ (^^.) + I^I^ ^ ^Rdiag ( ( 5^\^^ P/^R (A 



ypV _ g. /"-^ V^ ' ^ J_ 7qV_g. V" ' " y/_ 7q'p'-_lii /pV _ It C37^ 



Transforming the last result to the annihilation-creation operator correlation matrix we 
indeed obtain a well known formula for the occupation number of a harmonic oscillator 
in the equilibrium (thermal) state at inverse temperature j8 

Zf^ =i(^A,Z^J + f^+iZP}'-iZ^'j (38) 
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This solution is independent of the coupling operators Zl,Xr and detailed structure 
of the bath correlation functions, however, they should satisfy the KMS condition. 
Therefore, the same correlation matrix is obtained for a general harmonic oscillator 
chain with heat baths of equal temperatures. We prove that by expanding the Lyapunov 
equation (35) 

Mi'q' + (Mi'^y = Mf P'ZP'P' + i^(Zi'py +ZP'p'(M|''Py + ^Z^'P'^, (39) 

= Mfp'zi'p' + iazi'i' - izp'p', 
o = zi'p' + (zq'p')T_ 

These equations must be solved consistently with the symmetry conditions (Z'l'i )^ = 
Zi 1 and (ZP P )T = ZP p . From the last equation (39) follows the antisymmetry of the 
matrix Z^ p = Z^ p + ^1. Further, by simplifying the bath matrices for the case of equal 
temperatures 



M'l'i' = f {x^®x^ +x^ ®x^) diag Ul + eP^A T^ (Xj) 
M^'P' ='f{x^^x^+x^^ x^) diag ( ( ^^^^ ) r^R {Xj 

it is straightforward to see that correlations (37) solve the equation (39). 
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Figure 1. Comparison of numerically calculated thermal conductance and analytic solutions in the high 
temperature (47), and the low temperature (46) limits and the small coupling limit (54) for a homogeneous 
chain. We remark that analytic solutions for high and low temperature limits are asymptotically correct for 
all coupling strengths, which can be seen from the middle graph (e' == 5 x 10^^ and ohmic bath), whereas 
the small coupling solution is valid for all correlation functions, as shown on the right graph (e' = 10^^, 
sub-ohmic bath). The left graph is plotted for e' = 10~^ and ohmic bath. Other model parameters are: 
«==20, (7L-rR)/(rL + rR)=0.01, coo = 10, «c = 2/3, e[ = e^ = e. 
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3.2. Homogeneous chain - general solution 

In this section we provide a general solution for the homogeneous oscillator chain 
with ohmic bath correlation function F^ (ft)) = (fi^i-i - ^^ the case of fixed boundary 
conditions and equal, symmetric coupling Cl = Cr = £ it is possible to find an exact 
explicit expression for the correlation matrix. Straightforward but tedious calculation, 
in which we exploit the particular form of the bath matrices, shows that the solution 
of equations (39), rewritten in real space coordinates, can be sought in the form of the 
ansatz 

= sign(/ - j)z\i-j\ , Zn = zo = 0, (41) 

Z\i+j-2\-az\i+j-i\+Z\i+j\, ifi + j<n + l 

^fj = [Udiag(r^)U]y + e <! -Z\2n-i-j\ +aZ\2n-i-j-l\ - Z\2n-i-j-2\, if i + j >n+l 

0, otherwise 

We defined the bath functions 

' 2 I exp(/SLAj) - 1 exp(j3RAy)-l I 

and the constant a = 2 + (o^/(ol. The first part of the correlation matrix ZPP in equation 
(42) has the same form as in the equilibrium case, except the temperatures can be 
different. On the other hand, the second part is calculated from the equations (39) and 
the ansatz for Z^^. The Lyapunov equation is now simplified to a linear system of n — 1 
equations for the unknown coefficients zi,...z„_i. For a large number of oscillators 
n — )■ oo the covariance matrix is connected to a solution of the continuous fraction 
equation 



a + \/a^ - 4 (o} ,,^, 
^ , a = ^+a. (43) 



The heat current J = CO^zi is then given by 

1 , , 2(;c0i-0i) 1 "^^ 

Zi = ^2x^,+^2 + h + ---^n-l)= "^2/ +^2^L<^k. (44) 



k=l 



where the values ^1 , ^2, • • • ^n are calculated as follows 



" 2£f 7 . f nj \ . f knj \ 

^k=2^ —-J sm —— sm —— . (45) 

p'^n+l \n+lj \n+lj 

In the high-temperature limit the expression for the heat current is simplified to 

ft)2 

Aigh^^(7L-rR). (46) 
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This is an expected behavior since for high temperatures the thermal conductance 
Gih = J/{Ti^ — Tr) should be independent of the temperature, which follows from clas- 
sical mechanical considerations. In the low temperature limit we expect that thermal 
conductance will go to zero. We can simplify the thermal conductance for large on-site 
potential (coq ^ cOc) and small relative temperature difference \Tr — Ti^\/{Ti^ + 7r) <^ 1, 
to 

20)2 0)2 g-2a^,/(?L+rR) 

ex {Tl + Tr^ 



3.3. Recursive solution 

For a general spectral function of the form r^(tt)) = sign(ft)) ^^ (B(o)-i ^^'^ equal 
couplings it is possible to find a recursive solution for the correlation matrix in terms 
of a perturbative ansatz 

CO CO CO 

2P'p' = y £2J2P'p'(2;)^ 2i'i' = y £^JZ'^''^'^^J\ ZP'P' = y e2j+i2P'p'(2j+i)_ ^43^ 
j=o j=o 7=0 

The form of the ansatz can be guessed if we rewrite the equations (39) order by order. 
Starting from the lowest order we have: 

4t^'^ = m^ ((1 +e^-^^) T^i^j) + (1 +efe^^) rfe(A,)) = ?, (49) 

4''^'^ = ^ ((1 +e/^^'.) rH^j) + (1 +e/^^'.) rfe(A,)) = §., 

^j,k -(i-i-ij ') xj-xi ■ 

Each subsequent order is then calculated from the previous one using the recursion 
relations 

4''^"^ = Xihj { [M'i'P'zqV(-i)] . , _ [Mi'p'zqV(-i)],^^.} , (50) 

2P'p'(n) ^ 2Mi'p'Zi'p'(""^) +QZi'i'(""^), 
Z2'P'("+1) = i^^ {[M'i'p'zP'p'(n)]^.,+ [M'i'p'zPV(n)],^^.} . 

We can show that in each order for arbitrary initial condition Z^p^"^^) respecting the 
antisymmetry condition one can calculate the next order of of the correlations ZPP^") and 
Zii("), which also satisfy the symmetry conditions and from which the next correction 
for the coordinate-momentum part of the correlation matrix can be calculated uniquely. 
There might be problems on the diagonal due to the difference of the eigenvalues, but the 
special form of the bath matrix M^P , which is nonzero only when the sum j + kis odd, 

ensures that all diagonal corrections apart from ZPP^*^) and Z'^'^^^' vanish. We emphasize 
that the solution (50) is usefuU only for Cl = £r = £ < 1- 
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In the first order in tlie coupling e it is possible to find simple explicit results for the 
kinetic energy profile and the heat current. We connect the local kinetic energy to the 
local temperature at the same site 

7}- = ZP;P, (51) 

which is in principle justified only if the system obeys the condition of local thermal 
equilibrium. This is definitely not the case for harmonic systems, so our "temperature" 
field Tj which is simply defined by (51) should be understood only as a technical term. 
We find constant temperature profile for long chains and large on site potential (Dq^ (O 

f+2 1 /wo(l + e^La)oj (Oo(^l+e^R^^ 
4 ^41 exp(j8La^) - 1 exp(/3RC0D) - 1 



1-,'r,' 



^11 - ^ - 7 1 ^.^o „. ^ + „„.'; o „. ^ I • (52) 



We get the same form of the matrix X^ p as in the ohmic bath case, hence, the heat 
current in the weak coupling limit is 



n 



2 J^ \ n iv-1 



In the large on-site potential limit (tt)o ^ ft)c) the current is simplified to 

2l^"/2l i„/2-^ 4 (^exp(j8L«o) - 1 exp(/3RCOo)-lj- ^'^^ 

Comparison of analytical solutions in the low-temperature (46) and high-temperature 
(47) limits and small-coupling limit (54) with the numerical results for the thermal 
conductance is shown in Fig. (1). 



3.4. Disordered chain 

The homogeneous harmonic chains exhibit ballistic transport. Either one of the two 
different effects is needed for a sub-ballistic transport, say of a finite (non-zero) tem- 
perature gradient and finite (non-infinite) heat conductivity in the thermodynamic limit, 
namely: the interaction, which causes the dissipation of the phonons, or the disorder, 
which localizes the phonons. As interaction represents a whole new spectrum of prob- 
lems and phenomena (see e.g. [22]), we shall focus in this subsection to disorder in 
harmonic chains. However disorder often leads even to a sub-diffusive (insulating) be- 
havior in the thermodynamic limit. Disorder can be introduced by (i) randomly changing 
the masses, (ii) the on-site frequencies, or (iii) the coupling constants. We shall confine 
ourselves to the first type (i), sometimes referred to as isotopic disorder. The masses are 
chosen randomly nij = mo(l + S^j), where niQ denotes the mean mass, 5 E [0, 1] denotes 
the relative width of the deviation, and ^j is a random number uniformly distributed in 
the interval [-1,1]. In the literature both, the isotopic disorder [6] and the disorder in the 
coupling constants [17], have been studied in the context of quantum Langevin equation. 
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Figure 2. Averaged heat current scaling for different bath spectral functions (ohmic, Sub-ohmic, Super- 
ohmic) and boundary conditions (k' = 0, 1). The first panel (top) is for Oo = 10, the middle panel for 
(Ho = 1, and the bottom panel for Oo ~ 10^^ (almost gapless). The thermodynamic behavior for nonzero 
(0 does not depend on the properties of the baths and varying model parameters {coo,k'). The estimated 
asymptotic behavior is ^ n^'"*. The statistical error is approximately the size of the symbols. Model 
parameters: 5 = 0.8, © = 2/3, e[ = 0.3, eR = 0.1, j3l = 0.1, j3r = 1.0. 
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Figure 3. Averaged temperature profile (51) for different temperatures and on site potential and bound- 
ary conditions. Average the is made over 1000 realizations. Model parameters: 5 ~ 0.3, CO = 2/3, e^ == 
0.3, e^ = 0.1, j3r = 1.0 and ohmic bath spectral function. 



A finite temperature gradient and thermal conductivity in the thermodynamic limit have 
been observed [17]. 

Here, by employing our formalism, we show some new numerical results which differ 
from previous numerical conjectures [16, 6, 17]. Investigating the disordered harmonic 
oscillator Hamiltonian with free boundary conditions by using the quantum Langevin 
equation has been found that the thermodynamic behavior of the system depends on the 
properties of the heat baths [16]. However, using our Redfield operator space formalism 
we find heat current scaling independent of the bath spectral functions and coupling 
parameters. Numerical calculation of the heat current for different bath correlation 
functions, coupling strengths and boundary conditions is presented in Fig. (2). Note that 
in order to obtain reliable data the heat current and the temperature profile are calculated 
as an average over many realizations of disorder. For the gapped cases, COq 7^ 0, we find 
asymptotically, as n — > 0°, anomalous sub-diffusive power law scaling 



J(n) <x n 



where a = 1.4. 



(55) 



On the other hand, data seem less conclusive for the gapless case (in order to avoid 
numerical instabilities with our method we have chosen G)q = 10^^) where the scaling 
J{n) is either sensitive to the bath correlation function or the accessible system sizes n 
are still too small, and the behaviors J{n) are not yet asymptotic. In Fig. (3) we plot the 
temperature profiles for the ohmic baths and different on-site petentials and boundary 
conditions. It is perhaps remarkable to observe non-monotonicities in the profiles which 
are possible due to absence of local thermal equilibrium. 
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4. SUMMARY 

We generalized the method of the bosonic quantization in the Fock space of operators 
(third quantization) to a general markovian master equation. Lyapunov equation for 
the coordinate-momentum correlation function (the covariance matrix), which fully 
determines the non-equilibrium stationary state, has been obtained using the concept 
of left and right multiplication maps in the operator space. 

The formalism that we had derived was used to investigate the heat transport in the 
one-dimensional quantum harmonic oscillator chain with modified Redfield dissipators. 
An additional approximation turned out to be necessary in order to obtain consistent 
heat baths, which yield - in the case of equal temperatures of the reservoirs - the correct 
Gibbs state as the steady state solution of the master equation. For the homogeneous har- 
monic chain two analytical results for the covariance matrix were derived. The first one, 
given by equations (41), is valid for the ohmic bath and equal, possibly large, strengths 
of the coupling. The second is given in a perturbative, recursive form (50) and is valid 
for general spectral function and again equal couplings. This solution involves com- 
plicated expressions for higher orders, therefore it is useful only in the small coupling 
limit, where already the first correction is adequate. Simplified expressions for the heat 
conductance and the temperature profile were derived. As expected from the classical re- 
sults and previous quantum calculations [6] for the ordered chain we obtained vanishing 
temperature gradient and constant heat current in the thermodynamic limit. In the more 
interesting disordered case we find heat current scaling independent of the coupling, the 
bath properties, the on-site potential (as long as the phonon spectrum is gapped), and of 
the boundary conditions. The asymptotical decrease is faster than \/n (Fourier law), in 
fact it seems asymptotically oc n^^-^, which corresponds to an insulator. 
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